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Abstract. Integral representations are considered of solutions of the inhomo- 
geneous Airy differential equation w" — zw = ±l/ir. The solutions of these 
equations are also known as Scorer functions. Certain functional relations for 
these functions are used to confine the discussion to one function and to a 
certain sector in the complex plane. By using steepest descent methods from 
asymptotics, the standard integral representations of the Scorer functions are 
modified in order to obtain non-oscillating integrals for complex values of z. 
In this way stable representations for numerical evaluations of the functions 
are obtained. The methods are illustrated with numerical results. 



1. Introduction 
Airy functions are solutions of the differential equation 

Two linearly independent solutions that are real for real values of z are denoted 
by Ai(z) and Bi(z). They have the integral representations 

Ai(z) = i J °° cos (zt + ±t 3 ) dt, 

(1.2) 

Bi(z) = ~ / °° sin (zt + i*3) dt + i J Q °° e^ t3 dt 

where we assume that z is real. 

In this paper we concentrate on so-called Scorer functions ((?], [l^]), which are 
particular solutions of the non-homogeneous Airy differential equation. We have 

(1.3) 

w" — zw = —1/tt, with solution Gi(z) = — / sin zt + —t 3 I dt, z € K, 



3 



and 

1 f 00 _1 3 

(1.4) w"-zw = l/ir, with solution Hi(z) = - / e zt 3 * dt, ze 

n Jo 

Initial values are 
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Gi(0) =|Hi(0)= ±Bi(0) = ^Ai(0) = 



(1.5) 



^ W7W(fj 



Gi'(O) = iffi'(O) = ±Bi'(0) = -^Ai'(O) 



From (1.2), (O) and (O) it follows that 



(1.6) Gi(z)+Hi(^) = Bi(z). 

In the next section we give contour integrals from which representations of 
Ai(z), Bi(z) and Gi(z) follow for complex values of z. Just like Ai(^) and Bi(z), the 
Scorer functions Gi(z) and Hi(z) are entire functions. 

A survey on computational aspects of special functions, including information 
on Airy functions, can be found in |l(| has a public web site that includes an 
extensive treatment of Scorer functions. For complex values of z the Airy functions 
are available in the Bessel function algorithms of |2j ; see also pj and |l4| . Computer 
algebra systems as Maple and Mathematica also have Airy functions available. The 
Scorer functions are considered in where coefficients of Chebyshev expansions 
are given for real z. Asymptotic expansions for Hi(z) are given in [Q and The 
paper by Scorer (|ll|) gives tables to 7 decimals of Hi(— z) and Gi(z) for < z < 10. 

Efficient algorithms for computing the Scorer functions in restricted domains of 
the complex plane can be based on Maclaurin series and asymptotic series. These 
domains can be bridged by using the differential equations or the integral represen- 
tations. 

In [^| methods were presented based on the differential equations, which are set 
up as boundary-value methods, providing stable algorithms for all values of z. 
The purpose of this paper is to give stable integral representations for Gi(z) 



and Hi(z). We modify the integrals in (|l.3| ) an d (IT) by using methods from 
asymptotics. As in [[|, the resulting integrals can be used for any value of the 
complex parameter z. We also indicate how similar methods can be used for the 
Airy function itself. 

2. Asymptotic properties of the Airy and Scorer functions 

We need a few properties of the Airy and Scorer functions. More information 
can be found in jl], pC] , [pr| and [jl7). In particular Q, Chapter 11, discusses 
numerically satisfactory solutions of the differential equations for Gi(z) and Hi(z). 
The asymptotic properties of the Airy and Scorer functions are important in this 
discussion. 

We write, as in [pr| , 

(2.1) M (z) = Ai(z), Aii(z) = Ai (V 27 "/ 3 z) , Ai_i(z) = Ai (e 27 "/ 3 z) . 
We have the representations 

2j7ri/3 p 

(2.2) Ai j (z) = ^— e-* t+ 3* dt, i = 0,±l, 

where the contours Cj are given in Figure 2.1. Because J CoUClUC e ~ zt +k t3 dt = 0, 
we have the following linear combination of three solutions of ([Ll]) : 

(2.3) Ai(j?) + e- 27ri / 3 Ai!(z) + e 2 " /3 AU(z) = 0. 
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Figure 2.1. Three contours Cj of integration for the Airy integrals in ( |2.2| ), and 
sectors Sj where Aij(z) are recessive. 



The first integral in (1.2) follows from deforming the contour Co in (2.2) into the 
imaginary axis. The function Bi(z) can be written as 

(2.4) Bi(z) = e" /6 Ai_i(z) + e'^Ai^z) 

and the second representation in ( |1.2| ) follows by deforming contour C\ into the 
positive imaginary axis and (— oo, 0], and the contour C_i into (— oo,0] and the 
negative imaginary axis. 

The function Aij(z) is recessive at infinity in the sector Sj,j — 0, ±1, the function 
being exponentially small at infinity along any ray interior to this sector. On the 
other hand, Aij(z) is dominant at infinity in the sectors Sj-i and S^+i (the suffix 
j is enumerated modulo 3), and is exponentially large at infinity along the rays 
interior to these sectors. Bi(z) is dominant at infinity within all three sectors Sj. 
A pair of Airy functions comprises a numerically satisfactory pair at infinity within 
a sector if only one function is dominant. For example, the pair {Ai(z), Bi(z)} 
comprises such a pair only in So (and on the negative real axis, where none of the 
two is dominant, but where the phases in their oscillations differ by ^ tt). 

If one of the functions Gi(z), Hi(z) is computed, we may use (1.6) to compute 
the other one (we assume in this paper that Bi(z) and all other Airy functions are 
available), but we need to know if (1.6) is numerically stable for that computation. 
For example, because Gi(z) is only of algebraic growth in So (as we will see soon), 
we cannot compute Gi(z) from ( fL6| ) in So- 

For the Scorer functions we have the following asymptotic expansions (cf. |jlT|] , 
431-432): 



(2.5) 



Hi(z) ~ - — 

7TZ 



(3a + 2)! 
s!(3z 3 ) s 



oo, 



|pM-*)I < 



(2.6) 



Gi(z) 



1 

7TZ 



(3s + 2)! 



z 3 ^ s!(3z 3 ) s 

s=0 y ' 



OO, |ph Z\ < -7T 



5, 



5 being an arbitrary positive constant. For (|]^) the domain for ph z is not given 
in [nTl, but it follows from the same method mentioned for Hi(z) in pT|, p. 432. 
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In other parts of the complex plane we cannot represent the Scorer functions by a 
single expansion with leading term 0{X/z). 

From the results in (2.5) and ( |2.6| ) and the dominant asymptotic behavior of 
Bi(z) in all sectors Sj, j = 0, ±1, we conclude that (cf. ( |l.6| )) Bi(z) is a dominant 
term for Hi(z) in So and for Gi(z) in Si U S-%. It follows that we need algorithms 
for the computation of Gi(z) for z £ So and for Hi(z) for z £ Si U S-i (where 
the asymptotic expansions fl2.5| ) and (^) are valid). The relation in ( |l.6| ) can be 
used for computing the functions in the complements of these domains (where the 
functions have the dominant behaviour of Bi(z)). 

A further reduction of domains follows from the connection formula [] 



(2.7) 



Hi(z) = e ±27 "/ 3 Hi (ze ±2 ™/ 3 ) + 2e^ 6 Ai (ze^/3 



To prove this relation observe that the first term in the right-hand side satisfies 
the differential equation for Hi(z), and that, hence, that term can be written as a 
linear combination of Hi(z) and solutions of the homogeneous equation; the initial 
values in (1.5) can be used to i dentify these solutions. 

For example, we can use (2.7) with the upper signs for z in the sector \tt < ph z < 



|7T. Then, Hi(z) can be expressed in terms of Hi in the sector —it < ph z < — |7T 
plus an Airy function in the sector — 7r/3 < ph z < 0. We see that both functions in 
the right-hand side of (|]7]) are not dominant in the respective sectors, and, hence, 
this representation is stable. 

A similar connection formula for Gi(z) reads 

Gi(z) = e ±2 ™/ 3 Gi (ze ±2 ™/ 3 ) + e^^Ai (ze^W 3 

This formula is of no use in the sector So because both functions in the right-hand 
side are dominant, whereas Gi(z) is of algebraic g row th at infinit y w ithin So- A 
better formula for z £ So follows from combining (1.6), (2.4) and (2.7) (the latter 
twice, with upper and lower signs). This gives 



(2.8) 



Gi(z) 



e 2 " /3 Hi 



(ze 2 ™/ 3 ) 



e -2W3 Hi 



(ze- 2 ™/ 3 )" 



For z £ So the arguments of the Hi— functions are in S±i, where these functions 
have expansions that follow from (2.5). 
Because of 



(2.9) 



Hi(x — iy) = Ki(x + iy) and Gi(a; — iy) = Gi(x + iy) 



we can concentrate on non-negative values of the imaginary part y of the argument 
z = x + iy. 

Conclusion. The principal domain of interest for the Scorer functions is the 
sector |7T < ph z < tt, wher e w e conc entrat e on Hi(z). For z in other sectors, 
and for Gi(z), the relations (1.6), (2.7) and (2.£) are numerically stable for the 
particular cases. 

For a summary of the results of this section and algorithms we refer to Section 4. 



1 With thanks to the referee. 
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3. The construction of non-oscillating integrals 



We modify the integral in (1.4) such that stable algorithms can be based on the 
new integral for z in the sector |7r < ph z < it. 
We write 

(3.1) <fi(t) = -t 3 — zt, t = u + iv, z = x + iy. 

Then the real and imaginary parts of cf)(t) — r (u, v) 4- i(fii(u, v) are given by 



(3.2) 



(f>i{u,v) 



±u 3 

U V - 



UV 
3 v 



XU ■ 
XV - 



■ yv, 

yu. 



We are interested in a path in the complex t— plane on which (j)i(u,v) is a con- 
stant, and the path should start at the origin, as the integral in (1.4). Such a path 
is defined by the equation 

(3.3) u 2 v — -v 3 — xv — yu = 0, 



and we need real solutions of this equation. 

As summarized at the end of Section 2, we need to consider z— values satisfying 
27r/3 < ph z < 7r. If z < we can integrate along the positive t— axis. For other 
values of ph z in the present range the path of integration will be deformed into a 
curve C that is defined by ( [3.3] ) and that runs from the origin to +oo; see Figure 
3.1. 




Figure 3.1. The contour C is defined by (3.4). T he le ft figure is for z inside the 
sector 27r/3 < ph z < tt, and the cur ve is defined by (3.4). If ph z — 2ir/3 the curve 
in the right figure is defined by (3 



We solve the cubic equation (3.3) for v. The solution that passes through the 
origin can be written in the form 



(3.4) 

where 

(3.5) 



v(u) = 2 



1 



sin ■ 



P r 1 i 

arcsm — == g 0. — 7T 
2 



P 



:yu, Q 



x, 



with u > 0, x < 0, 3x 2 > y 2 . 
To show this, we introduce t 



At 3 



v/(2y/Q). Then ( |3.3|) can be written in the form 
P 



•3t= — 



<Q 3 



It is not difficult to verify that, if u > 0,x < 0,3a; 2 > y 2 , then the modulus of 
the right-hand side is not larger than unity. Replacing the left-hand side with 
— sin(3 arcsin(i)) gives the solution in (|3~5|). 
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It follows that 
(3.6) 

where (f> r (u,v) is given in ( |3.2| ), v(u) in (3.4) and 

d[u + iv(u)] 



1 f°° 

m(z) = - / e-^ v ^h{u)du, 
n Jo 



(3.7) 



du 



= 1 



dv(u) 



2uv — y 



If ph z = 2n/3 then y = — x^/S, andx < 0, equation (3.3) can be solved explicitly. 
The two solutions are 

(3.8) v = uV3 and u — . 

In this case the path of integration C runs from the origin to the point 



t = uq + iv = \/—x/2 + i\J '-3x/2 

along the line v — uv3, and for w> uo the path C follows the hyperbola defined 
by the second solution given in (|3.8|); see the right figure in Figure 3.1. The point 
to = uq + ivp = \fz = \/ x — ixy/3, x < 0, is a saddle point of the function <f>(t) 
defined in Q ). 



(3.9) 



We can also solve (|3.3| ) for u, which gives 
y-R 



R 



1 



y 2 + Av 2 (x H — u 2 ) 



2u V 3 

where the square root is non-negative. This solution should be used for < v < 
Vi, Q < u < U\, where 



vi = vi(x,y) = J ^ (-x -^x 2 - y 2 / 3 )" 



the smallest positive u— value for which R — and iti = y/(2vi). For < u < 
ui, it > iti, we use u = (y + R)/(2v); see Figure 3.1. 

When integrating with respect to v, the integral in (|3.6]) can be written as 



(3.10) 



Hi(z) = i 



,-4> r (u (v),v) 



h(v) dv 



h{v) dv 



where w) is given in (B^h, tx ± (v) = (y ± #)/(2«) (cf. (Bj)) and 



(3.11) 



d(u + du 



dv 



dv 



2uv — y 



Remark 3.1. For the sector < ph z < ^ir we can use a similar method, although 
we don't need to consider this sector. If y 2 ~ 3x 2 > the quantity R of (|3 . 9f) is 



defined for all values of v, and the first equation in (3.9) defines a path C going 
from the origin to oo exp(27ri/3). Be caus e there is a hill at oo cxp(7ri/3), we can- 
not replace the integration path in (1.4) by C. We need an extra integral from 
oo exp(27rz/3) to +00, and that integral gives an Airy function; see Figure 3.2. In 
this way we obtain 



(3.12) 



Hi( z ) = ~ J e-^ u ' v) h{v) dv + 2e- 7 "/ 6 Ai (2 



a -27ri/3 
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with h(v) given in (3.11) and the relation between u and v given in (3.9). We see 
that the first term in the right-hand side of (2.7) with upper signs corresponds with 
the integral in ( [3.12] ). 




Figure 3.2. The two contours for the integrals described in Remark 3.1. 



Remark 3.2. Comparing the representations of Hi(z) in ( p.6|) and (3.12) we observe 
that the Airy function in ( 3.12| ) disappears as z crosses the half line ph z = 27r/3. 



On that line the argument of the Airy function Ai [ze 2 ™/3) becomes positive (see 
(3.E)). Thus the dominance of Hi(z) over Ai (ze -27 ™/ 3 ) is maximal at this line, 



this is therefore the place where the exponentially-small contribution is "switched 
on" . This aspect is connected with the Stokes phenomenon in asymptotics, and the 
half-line ph z — 2n/3 is a Stokes line for Hi(z); see [ p"2| . 

3.1. The Scorer function Gi(z). It is convenient to have a direct method for 
Gi(z) that is not based on results for the Hi— function, especially if z is positive. 
We give only a few details on this case. The method can be used for the sector 
< ph z < 2tt/3. 

The first step is to replace the sine function by two exponentials. Then we obtain 

1 f°° 
(3.13) Gi(z) = [Gi+(z) - Gr(«)] , Gi^z) = / dt, 

where ifj(t) = zt + \t 3 = Vv(^> v ) + ^iC^i w )> with 



(3.14) ' " • 



ip r (u, v) — — uv 2 + xu — yv, 
ipi(u, v) — u 2 v — ^v 3 + xv + yu. 



,2„, 1,,3 

— r 

The path of integration for Gi ± (z) is found by solving the equation tp r (u, v) = 0. 
For Gi + (z) the path runs into the valley at oo exp(7ri/6). The result is 

roo , ..1 ,.2 , _ 

(3.15) Gi+(z)= / e~^ v) g(u) du, g{u) = 1 + i— = 1 



du 2uv + y 

For Gi~ (z) the integral along [0, oo) can be replaced to a path along the half-line 
with ph t = — 7r/6. In this valley no real solution of ip r {u, v) = is available, and 
we take a path that first runs into the valley at oo exp(— 57ri/6) and then returns 
to the valley at oo exp(— ni/6). This introduces an Airy function, and we obtain 

/•OO 

(3.16) GP (z) = e -Mu,v) g ^ du + 27r Ai(z), 
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where g(u) is as in ( |3.15| ). Adding the results in (3.15) and (3.16), we obtain a 
simple non-oscillating integral plus an Airy function: 

(3.17) Gi(z) = — / e -^ {u ' v) g(u) du + iM{z), 



where the relation between u and v is given by 



-y + \Jy 2 + 4u 2 (x + \u 2 ) 
2m ' 



If z is real and non-negative (3.17) becomes real. The term with the Airy function 
is canceled by the imaginary contribution of g(u)/i. The remaining integral in (3.17) 
should be modified in this case. The contour runs from the origin to the saddle 
point iyfx, and from this point into the valley at oo exp(7ri/6). Integrating with 
respect to v, we obtain the real representation 



(3.18) 



Gi(x) = - 

7T 



cfo + / e*" -*" 3 dv 



x > 0. 



4. Numerical illustrations 



We give some numerical results which serve as demonstration of our method. 

In order to evaluate the Scorer functions in the whole complex plane, we need to 
select software for the evaluation of the Airy functions of complex arguments and 
for the quadrature of real functions over an infinite integral. For the first purpose, 
we use the public domain subroutines ZAIRY and ZBIRY by D.E. Amos [|| and 
for the semi-infinite integral we use the automatic adaptative integrator DQAGIE 
by R. Piessens. All these codes can be retrieved from the SLATEC public domain 
library Q (see also GAMS: guide to available mathematical software ||). 

The connection formulae given in Section 2 together with the non-oscillating 
integrals derived in Section 3 can be used to evaluate Hi(z) and Gi(z) in the whole 
complex plane. By using the integral representations for Hi(z) in the domain 7r < 
ph z < 2tt/2> the following stable algorithm can be considered: 



Algorithm 4.1. (via (S.t j) 



Whenever Im (z) < 0, use eg. pLa). 

(2) 

If z £ S{ then obtain Hi(z) by quadrature. 



Ifze S[ 1} U S n obtain Hi{z) via (IO 



Obtain Gi(z) everywhere in the complex plane by using (2.8) 



where is the sector 7r/3 < phz < 27r/3 and s[ 2 ^ is the sector 27r/3 < ph < it. 

However, in the fourth step two integrals for Hi will be needed. Thus, the 
following stable scheme is expected to be more efficient provided fast algorithms to 
compute the Airy functions Ai and Bi are available: 
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Algorithm 4.2. (via $Ufy and (p.lty ) 



Whenever Im (z) < 0, use eq. \2.!\ ). 

( 2) 

If Z € S\ then obtain Hi(z) by quadrature. 
If z € Si \JSq then obtain Gi{z) by quadrature via ( 3.11 ) or ( 3.18 ). 
For the remaining cases, apply eq. IU_ 



The second algorithm is preferred in most circumstances. However, we have 
experienced that the computation of Gi(z) when Re(z) > and Im(z) — ► is 
more efficient when the first algorithm is considered. Probably, the best numerical 
strategy is a combination of both algorithms, together with the use of asymptotic 
expansions for large \z\ and series expansions for small \z\. The best strategy may 
also depend on the choice of the quadrature rule. This numerical study lies beyond 
the scope of the present paper. 

We end this section by showing numerical results. Of interest are the sectors 
where (2.5) and (2.6) are valid, namely, S^ 2 ' for Hi(z) and IJSo f° r Gi(z). Of 



(2) 

particular interest are the regions So for Gi(z) and S^' for Hi(z) since in this case 
we compute the functions directly by quadrature and the corresponding integrals 
can be compared with asymptotics if z is large. We are giving explicit results with 
an accuracy of 8 digits for Hi(z) which are compared with the asymptotic expantion 
( pH| ) up to order l/z w whenever this is possible. 

The results in Table 4.1 are obtained by means of a Fortran program coded 
in double precision arithmetic in which the integral ( |3.6| ) is evaluated. Together 
with the results, we show the number of integration steps needed to attain an 
accuracy of 8 digits. The results from asymptotics (shown inside brackets) are seen 
to coincide with those from the integral ( |3.6| ) for \z\ = 100, but for \z\ = 10 we 
observe discrepancies in the last digits which are due to the failure of the asymptotic 
expansion (z is not large enough). We also found agreement with the asymptotic 

(2) 

expansion for Hi(z) for large z in the sector S{ . For real negative x our results 
coincide with those given by Scorer jl5| . 

One sees that at the Stokes line ph z = 2tt/3 the quadrature requires more 
steps, as can be expected given the appearance of a discontinuity in dv/du at the 
maximum v. On the other hand, the faster convergence takes place when we are 
far from the Stokes line. Moreover, convergence tends to be slower as \z\ becomes 
smaller; as \z\ becomes larger the effect of the singularity in the derivative dv/du 
for ph z — |7r would appear at larger u, where the exponential in the integrand is 
smaller. 

Similarly, one can test the performance of the integral representation (plus an 
Airy function) for Gi(z) in the sector < ph z < 2tt/3, which should coincide with 
the results from the asymptotic expansion ( |2.6| ) for |ph z\ < 7r/3 and large z. We 
also find agreement with the results tabulated by Scorer for real positive x. One 



observes that the convergence of the integral representation (3.17) becomes slower 



as we approach the real line ph z — 0. With the quadrature DQAGIE we obtain a 



better performance using Algorithm 4.1 in this case 
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pn z = n 


pn z ~ 57T/6 


pn z = Z7T/3 


\z\ = 1 


Re (Hi) 


0.22066961 
{195} 


0.22331566 
{195} 


0.23477589 
{345} 


1 ~ 1 1 

\z\ = 1 


1m ) 


U 


d 11 OQA01 1 A-2 

{165} 


U.lot>Uo5y4 
{465} 


\z\ = 10 


Re (Hi) 


3.1768535 10"^ 
{75} 

(3.1768528 10" 2 ) 


2.7597145 10~ 2 
{75} 

(2.7597137 10~ 2 ) 


1.5948003 10~' z 
{135} 

(1.5947998 10~ 2 ) 


\z\ = 10 


Im (Hi) 





1.5859789 10~ 2 
{75} 

(1.5859786 10~ 2 ) 


2.7622751 10~ 2 
{225} 

(2.7622742 10" 2 ) 


\z\ = 100 


Re (Hi) 


3.1830925 10" 3 
{135} 

(3.1830925 10" 3 ) 


2.7566477 10" 3 
{165} 

(2.7566477 10" 3 ) 


1.5915526 10" 3 
{165} 

(1.5915526 10~ 3 ) 


\z\ = 100 


Im (Hi) 





1.5915439 10~ 3 
{165} 

(1.5915439 10" 3 ) 


2.7566500 10" 3 
{165} 

(2.7566500 10" 3 ) 



Table 4.1 The real and imaginary parts of Hi(z) in the sector tt < ph z < |7r. The 
result from the asymptotic expansion up to order 1/z 10 is shown inside brackets. 
The number of integration steps for each evaluation is shown within the braces. 

As a further illustration, we will give two plots for \z\ = 1 obtained by using Al- 
gorithm 4.2. The graphs show the real and imaginary parts of the Scorer functions. 
One observes the smooth connection between the different sectors in the complex 
plane. 




CD 




Figure 4.2 The real (solid) and imaginary (dashed) parts of Hi(z) (left) and Gi(z) 
(right) for \z\ = 1 and < ph z < tt. 
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5. Summary and concluding remarks 

As mentioned in the cited references (see JjJ and ||) the inhomogeneous Airy 
functions (Scorer functions) are used in several physical problems. We have used 
functional relations and derived integral representations of the Scorer functions 
that can be used for stable numerical computations for all complex values of the 
argument z. The integrals can be easily split up into real and imaginary parts. 

For the Scorer functions Maclaurin and asymptotic expansions are available. To 
bridge the gap between the areas where convergent or asymptotic series can be 
used, one can use the representations in this paper, although they can be used for 
all values of the argument. The tool one needs is a suitable quadrature method 
for computing real integrals on unbounded real intervals that converge very fast at 
infinity. We have illustrated the method by giving a few numerical results based on 
selecting a quadrature rule; we have not investigated an optimal choice of quadra- 
ture rule for computing the Scorer functions. 

We have shown how to handle oscillating integrals for a certain set of special 
functions, and this is quite instructive for applying the method to other functions. 
A similar method can be used for the Airy functions and another treatment can be 
found in ||l6f for modified Bessel functions of imaginary order. Still many special 
functions need to be considered in order to get reliable software, in particular for 
complex and/or large values of the parameters. 
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